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I. INTRODUCTION 


A. BACKGROUND 

In analyzing ocean acoustic pulse-propagation problems, 
accurate phase calculations must be performed in order to 
predict the acoustic signal at the receiver. Linear systems 
theory provides analytical expressions for analyzing the 
propagating acoustic field [Refs. 1- 4]. These well-known 
equations form the basis for the FORTRAN programs developed in 
this thesis. 

Program input is depth versus sound-speed data pairs where 
the speed of sound is an arbitrary function of depth. The 


ocean is viewed as a linear, time-variant, space-variant 


filter. The WKB approximation can specify this filter's 
transfer function [Refs. 1 - 4}. For the arbitrary sound- 
speed profile, no exact transfer function exists. The 


transfer function requires a method for evaluating signal 
phase at the receiver. 

Two phase evaluation methods are presented and contrasted. 
Each method calculates the phase of the acoustic signal fora 
specified horizontal range. The first method overlays a 
piecewise linear profile on the input sample values for the 
speed of sound at various depths. The theory of ray acoustics 
is used to calculate travel time and phase. The second method 


fits a smooth cubic spline curve to the input samples. This 


method solves the propagation problem using a system of three, 
first-order differential equations [Ref. 5]. Both methods 
make phase calculations using signal travel times. 
Additionally, each method is capable of identifying eigenrays 
or rays that directly connect the signal source to the 
receiver. 

Another method is introduced to validate results. Direct 
numerical integration is performed to calculate phase for a 
Single gradient, linear, sound-speed profile. A separate 
FORTRAN program implements this numerical integration routine 
providing a totally independent verification. 

An overview of the theory behind each method is presented 
in Chapter II. Computer simulation results are presented in 
Chapter III for various input sound-speed profiles. Tabular 
results of phase calculations are presented for each of the 
three analysis methods. Using input from the piecewise linear 
and the cubic spline/differential equation solution methods, 
ray traces of the propagating field are presented as a visual 


aid in interpreting the results. 


II. THEORETICAL BACKGROUND FOR THE 
EVALUATION OF PHASE INTEGRALS 
A. UNDERWATER ACOUSTIC PULSE PROPAGATION 

Ocean acoustic pulse-propagation models can be derived by 
using the principles of linear, time-variant, space-variant, 
systems theory and the physics of wave propagation in 
inhomogeneous media. Linear systems theory allows for the 
development of an ocean medium transfer function. An ocean 
medium transfer function that is based on the WKB 
approximation has been derived and is given by References 1 


and 2 as follows: 


~1{" k(¢)d¢ 
ECC ee ec JOO a sentyty-) (2.1) 


VTky Cy) | 


where 
fe is the frequency in Hz 


is the radial, spatial frequency in 
cycles per meter 


Yo is the source depth in meters 


a = 7 ¥1Ky (Yo) | (222) 


anf, 


ioGy) 1s the propagation vector component in 
the Y direction with units of radians 
per meter and is given by 


k, (y) 


° am{(f/e(y)]}° - £,°)" ,f£, < £/e(y) (2 #5 
F jan{f, - (£/c(y)1*)" ,f. > £/ewy) (2.4) 


c(y) is the sound-speed expressed as a function 
of depth with units of meters per second 


f, is the transmitted (input) spatial frequency 


in the Y direction at the source in cycles per 
meter and is given by 


( [ (£75) = eet eet ee (2.5) 
fy = 


- (f. - (fe) | Aeeeeee (2 a 
and 


Cy) is the speed of sound in meters per) Second at 


the source depth y,, that Toye; =se 2 


The plus (minus) sign in Eqs. (2.3) and (2.5) is chosen 


whenever y - yp > O (y —- y, <= 0). The manus (plus) sigieeen 


Eqs. (2.4) and (2.6) corresponds to the plus (minus) sign in 


Eqs. (2.3) and (2.5). This thesis deals only with propagating 


waves. Therefore, Eqs. (2.4) and (2.6) representing the 


generation of evanescent waves (1.e., decaying exponentials), 


will not be used. 


The evaluation of the phase integral in Eq. (2.1), namely, 


_y 
6 (y) = [rw (0) ad (22am 


has been attempted in past studies in which solution 
techniques included direct integration and binomial expansions 
fmef. 3). 

Computer programs using these techniques proved to be very 
expensive in terms of CPU time. This thesis compares four 
adifferent methods for evaluating the phase integral. The four 


methods are presented in the following sections. 


B. DIRECT INTEGRATION 
Computer code was developed to directly evaluate the phase 
integral given by (Eq. 2.7) using numerical integration 
routines from the International Math and Statistics Library 
(IMSL). Direct integration is a time-consuming technique that 
was applied only under the following constraints: 
- free-space acoustic propagation 


- linear sound-speed profile with a single gradient g 
in inverse seconds and 


- propagating signals that have not passed through a 
turning point in the medium. 


Expressions for the propagation vector component k,(y) of 
the phase integral are given by Eqs. (2.3) and (2.4). Fora 
linear sound-speed profile with a single gradient, the 
dependence of c(y) on depth y is given by 

CCY JM CU) t8Gi(yY =Yo) (2.8) 
where 


C(Y)) 1s the sound-speed at the source depth 
Yp meters 


g is the sound-speed gradient, and 


4 is the desired depth for sound-speed. 


Since the direct integration approach is limited to free- 
space, Single gradient, propagation problems, it is used only 
to validate the results of other solution techniques for 
several simple test cases. The results will be compared 
against the next two methods to be presented. These methods 
will apply the theory of ray acoustics [{Ref. 4 & 5] to find 
travel time. Phase calculations will be based on travel time 
calculations. 

C. PHASE CALCULATIONS BASED ON PIECEWISE LINEAR SOUND-SPEED 

PROFILES 

A FORTRAN program developed by Lim [Ref. 5] was the ray 
acoustics propagation code that was used to draw ray trace 
plots and to calculate travel time, in addition to the depth, 
angle of propagation and path length along a ray path as a 
function of horizontal range. This method applies ray 
acoustics to a piecewise linear sound-speed model of the ocean 
medium. Based on the input sound-speed versus depth data 
pairs, layers are defined in the ocean medium. Each layer has 
an upper and lower boundary at specific depths. A constant 
sound-speed gradient g is calculated for each layer. The 
sound-speed for any desired depth is computed using Eq. (2.8) 
with the appropriate gradient. 

In a constant gradient medium, ray acoustics theory allows 
calculation of travel time with closed-form equations. For an 
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incremental increase in horizontal range from the source, the 


angle of arrival for the propagating ray can be shown to be 


Baker. 9]. 

B(y) = cos” [cos@, - b*g(y) * rngstp] (2.9) 
where 

Bo is the ray launch angle, 

g(y) 


is the gradient in the layer at depth y, and 


rngstp is the incremental increase in horizontal 
range. 


The ray parameter b in Eq. (2.9) is given by [Refs. 4 & 5}. 


b = sinB, / C(Yp) C22 Os) 


Knowing the arrival angle B(y), the ray depth y is given by 
[Refs. 4 & 5] 


Cc 
Y = Yq (2.11) 





sinf (y) 3) 


0 
ao 
g | sing (y,) 


The arrival angle also allows calculation of travel time from 
(Refs. 4 & 5] 


T — 


tans (y) /2] 


= (22 ) 
S| tan[acy,) /2| 


Once the travel time (in seconds) is found for the desired 
horizontal range, total phase is given by 


Creme LT G21 3 ) 
The phase integral given by Eq. (2.7) represents a phase 
change in the depth (y) direction. The constant value of the 


propagation vector component in the radial direction allows 
for easy calculation of phase change in the radial direction. 
This value is given by 

Oe 27 fsa zeng (2.14) 
where hrzrng is the total horizontal range traveled. 
The phase in the depth direction is 

@ = 6, - &. (25) 
Phase expressed in radians 1S a modulo 2m function. The 
proper solution for the phase integral Eq. (2.15) must be 
expressed as a modulo 27 function: 


8, = modulo (8 


3 Mg ila, alr (2 Say 


y? 

This method (referred to as method 1 in Ref. 5) is capable 
of analyzing propagating rays at any horizontal range in 
either free space or a bounded medium. 
D. PHASE CALCULATIONS BASED ON AKIMA CUBIC SPLINES AND 

ORDINARY DIFFERENTIAL EQUATIONS 

This application of ray acoustics applies a more 
sophisticated treatment to the depth versus sound-speed data 
pairs sampled from the ocean medium. The data pairs are used 
to form a smooth Akima cubic spline sound-speed profile. 
Splines offer the advantage of using all data points in 
generating a profile, and they place no restrictions on 
spacing between data points. The Akima version of the cubic 
spline was chosen for its excellent ability to combat wiggles 


in the profile, that is, it suppresses oscillations that would 


Cause overshoots and undershoots in the sound-speed versus 
depth profile. 
This method uses ray acoustics theory to generate a system 


of three, first-order differential equations [Ref. 5]: 


Y, = Y> (2441/7) 
me ya) 
Yom — (2.18) 
bec’ (Y,) 
and 
1 
C) daa (2.19) 
be (y,) 
where 


Gives ehesderivative of the sound-speedswith 
Bese cemeOnacoul dtc Gepen y., 


Y; is the ray depth, 

Y> is the cotangent of the ray's arrival angle 
B(y), and 

Y3 is the travel time of the ray. 


Once solved, this system of differential equations allows 
phase to be calculated from travel time. The phase 
calculation is performed exactly as shown for the piecewise 
linear sound-speed profile. 

Like the previous method, the Akima Cube 
spline/differential equation solution can be used on free 
space and bounded media problems. 
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E. EIGENRAYS 
Eigenrays are propagating rays that exactly connect the 

sound source to the receiver. The FORTRAN propagation code 
developed can search for and identify eigenrays. The input 
required is 

- depth versus sound-speed data pairs, 

- source depth in meters, 

- receiver depth and range in meters, 


- angle step in degrees between possible eigenrays to 
be evaluated, and 


- allowed depth error idenavedsy =e Gelecer se 
Rays passing within the allowed error or tolerance y,,,,. of the 
receiver are identified as eigenrays. The eigenray mode can 
employ either the piecewise linear sound-speed profile or the 
Akima cubic spline/differential equation method for ray 
propagation. 

The eigenrays are found by trial-and-error. This search 
method was chosen after the IMSL DBVPMS program failed to 
solve the problen. The DBVPMS program is a differential 
equation solver that was applied to the system of differential 
equations given by Eqs. (2.17) through (2.19). The routine 
uses the shooting method to find eigenray solutions to 
boundary value problems. It was unable to converge to a 
solution. Solutions for the acoustic problems investigated 


are difficult for this algorithm because of 
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- long propagation distances, 


- the inhomogeneous ocean presents a _ continuously 
varying medium with discrete boundaries, and 


- the system of differential equations has non-constant 
coefficients whose values are a function of depth. 


When the eigenray mode is selected, rays that are not 
Betinaeted Oweretlected to within y.( Gf the recetver are 


discarded. 


alak 


III. COMPUTER SIMULATION RESULTS 


A. OVERVIEW 
The computer simulation results presented in this 
chapter perform phase calculations intended to: 
- show that ray acoustics can be used to _ evaluate 
the phase integral producing the same results as 
direct numerical integration, 
- show that the two methods of applying ray acoustics 
theory (presented in Chapter II) produce equivalent 


results within an acceptable tolerance, 


- compare speed versus accuracy for the two methods for 
a variety of test cases, and 


- demonstrate the ability to identify eligenrays using 
either of the ray acoustics methods. 
B. DIRECT INTEGRATION 
In the previous chapter, 8 (y) was shown as the solution 


to the phase integral in the depth direction. 


Re 
By(y) =] kK (o)ae (2am 


The straight-forward evaluation technique is the brute 
force approach using direct numerical integration. The phase 
integral will be evaluated using direct numerical integration 
for a Simple test case. The results will then be compared to 
phase values obtained using ray acoustics calculating phase 


values from travel times. 
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Conditions for the test case are: 
- source depth of 10 meters, 
- receiver depth of 100 meters, 


- linear sound-speed profile with a single, constant 
gradient of 0.016 sec and 


- four values of radial spatial frequency (FR values) 
will be evaluated. 


For each FR value, tabular results will show the 
corresponding ray launch angle B(y,).- The relation between 
the FR value and this launch angle is [Ref. 4]: 

B(yo) = sin’ [FR * c(y,)/f] (Goay 
where 

f is the frequency in HZ, and 

Giy,) 15 the speed of sound at the source depth. 

Table 3.1 shows the results produced by direct numerical 
integration. Tables 3.2 and 3.3 show the results for the 
Akima cubic spline/ordinary differential equation solver 
(ODE solver) and the piecewise linear approximation ray 
acoustics techniques, respectively. 

The parameter of interest is THETAY representing the phase 
change in the depth or Y direction. Table 3.1 shows these 0 
values to be approximately 92, 86, 75, and 56 radians for the 
four FR values chosen. Tables 3.2 and 3.3 contain four sets 
of entries, each corresponding to an FR value. These entries 
begin at 0 meters range and end at the range corresponding to 


100 meters depth; thus, solving the phase integral for the 0, 
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TABLE 3.1 
PHASE INTEGRAL SOLUTION BY DIRECT NUMERICAL INTEGRATION 


INPUT DATA FOR PHASE INTEGRAL EVALUATION 


Foe 250.0 HZ YL s 10.0 4 YU = 100.0 M RATIO = 0.9999 NFR # 5 
EABS = 0.000000 EREt = 0.000100 

YREF s 0.0 M CYREF = 1500.0 M/SEC G = 0.016000 (1/SEC) 

METMOD = 2 


EVALUATION BY IMSLIO ROUTINE DODAG W1TM IRULE = 2 


FR(CYCLES/M) BETAYO(DEG) RANGE AT ANGLE OF ARRIVAL TURNING POINT TURNING POINT TMETAR(RAD)  THETAY(RAD) 
yu BETAY( DEG) DEP TM(M} RANGE (M) 
0.0332944859 11.525 18.36009 11.536 O.37555E*06 0.45986E*06 0.38408E*01 «0. 822916 79386E +02 
0.0665889718 23,552 39.25213  . 23.576 O.14090E*06 0.21S10E*06 0.16423E+02 0. 86 S3AGO2GeE +02 
0.09988346576 36.826 67.62892 36.866 O.62692E*0S 0, 1SS22E+06 = 00.4 2324E+02 ~—s 0.75 378S7191E+02 
O.1521779625 52.049 119.80722 53.122 O.23574E*0S 0. 70527E405 0. 1 002SE*03 «= 0.5 6S 7 SK BGI IE +02 
TABLE 3.2 
PHASE INTEGRAL SOLUTION BY ODE SOLVER 

ERAYS = F PRIALL 251 

YO s 10.0 M DEPTH = 100.0 M  FREQC » 250.0 HZ CMAX = 1501.6 M/SEC 

FREOC/CMAX = 0.16669 CYCLES/M RATIO ® 0.9999 FRMAX = 0.16647 CYCLES/M NFRS © 5 

DLIFR = 0.033296 CYCLES/M RNGSTP » 1.0.M MRZRNG = 150.0 M VR ®= 60.0 M 


NOTE: FRMIN © DLTFR 


FR # Q.332946858861726E-01 CYCLES/M 


BETAO = 11.52457521170979 DEG 

RANGE (M) DEPTH(H) TRVLT (SEC) BETAY (DEG) BYEST THETAT(RAD) THE TAR (RAD) THE TAYCRAD) 
0.0 10.00 0.000000 11.525 O.13318E-03 0.0000000E+00 0.0000000E+00 0.0000000E*00 
18.4 100.00 0.061200 11.536 O.13318E-03 0.9613253E*02 0.3840846E*01 0.9229168E*02 


FR = 0.66588971763452E-01 CYCLES/M 


BETAO = 23.55170283980132 DEG 

PAHGE (CM) DEPTH(M) TRVLT(SEC) BETAY( DEG) BTEST THETATCRAD) THE TAR CRAD) TME TAY (RAD) 
0.0 10.00 0.000000 23.552 0.26636E-05 0.0000000E*00 0.0000000E*00 0.0000000E*00 
335 100.00 0.065420 23.576 0.26636E-053 0.1027611E+03 0.1642273E*02 0.8633840E*02 


FR © Q.998834576465178E-01 CYCLES/M 


BETAO « 56.82460911646773 DEG 

RANGE (M) DEP TH(M) TRVLT(SEC) BE TAY (DEG) BTEST THE TATC(RAD) THE TAR (RAD) THE TAY (RAD) 
0.0 10.00 0.000000 36.824 0.399S°E-03 0,0000000E*00 0.0000000E+*00 0.0000000E*00 
67.4 100.00 0.074932 36.866 0. 39953E-05 O.1177023E403 0.4232376E*02 0.7537857E*02 


FR = 0.133177942352690E*00 CYCLES/M 


BETAQ = 53.04928575821315 DEG 
RANGE (M) DEP TH(H) TRVLT(SEC) BETAY (DEG) BTEST THETAT (RAD) THE TAR (RAD) THE TAY (RAD) 
0.0 10.00 0.000000 53.049 0.5327 1E-03 0.0000000E+00 0.0000000E+00 0.0000000E+00 
119.8 100.00 0.099829 Ss lze 0.S2271E-03 0.1S68262E*03 0.100°0S525E*03 0.5657368E*02 
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TABLE 3.3 


PHASE INTEGRAL SOLUTION BY PIECEWISE LINEAR APPROXIMATION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1°L4 


NDATA = 4G NUMBER OF GRADIENTS = 3 

Gt 1) = 0O.16000E-01 1/75EC GC 2) = 0.16000E-01 1/S5EC 

Meee = O0.16000E-01 1/S5EC 

ERAYS = F PRTALL = T 

YO = 10.0 M DEPTH = 100.0 M FREOC e« 250.082 CMAX = 1501.6 M/SEC 

FREOC/CMAX = 0.16649 CYCLES/M RATIO = 0.9999 FRMAX #® 0.16647 CYCLES/M NFRS = 5 

Metem = 0.025294 CYCLEG/H RNGSTP = 1.046 HR2ZRNG 2 USP ON (0) aya! YR 2 40.0 ™M 

NOTE: FRMIN = DLTFR 

FR = 0.322946485881726E-01 CYCLES/M 

PETAO » PP. 5245752170980 DEG 

RANGE (14) DEPTH(M) TRVLTCSEC) BETAY(DEG) BEST THE TAT(CRAD) THETAR (RAD) THE TAY (RAD) 
0.0 10.00 0.000000 re Dela sise—0s 0.0000000E*00 0.0000000E+00 0.0000000E+00 
18.4 100.00 0.061200 DS 26 0. Ps 1 5-05 O.9612E53E*02 0.3840846E*01 Oeeece  6se<02 

PR = 0.665889717634S52E-01 CYCLES/M 

PETAO = 23.55170285980134 DEG 

RANGE (M) DEPTH(M) TRVLT (SEC) BETAY (DEG ) BHEST THETATCRAD) THETAR (RAD) THE TAY (RAD) 
0.0 10.00 0.000000 25255¢ 0.26656E-03 0.0000000E*00 0.0000000E*00 0.0000000E+*00 
vi) es 100.00 0.065420 23.576 O.26636E-03 0.102761 1E+03 0.1642273E+02 0.8633840E*02 

FR = 0.99882457645178E-01 CYCLES/M 

BETAO = 36 .82440911644776 DEG 

RANGE(M) PEP THICHp TRVCI See) BETAY (DEG ) BTEST THETATCRAD) THETAR (RAD) THETAY (RAD) 
0.0 10.00 0.000000 36.824 Ono 72> 5E705 0.0000000E*00 0.0000000E+00 0.0000000E+00 
67.4 100.00 0.074932 36.866 Oo 7 Jo5e— 05 Col 77023E*03 0.46232374E*02 0. 7537257E*0. 

Beet Ostool7796252690E*00 CYCLES/M 

BETAO = 53.04928575821320 DEG 

RANGE (1) DEPTH(M) TRVLT (SEC) BETAY (DEG) BIEST THETATCRAD) THE TAR (RAD ) THETAY (RAD) 
0.0 10.00 0.000000 53.049 D.5s271E-03 0.0000000E*00 0.0000000E+00 0.0000000E*00 

119.8 100.00 0.099839 eee Une. ! LE=05 0.1568°6°0E*03 
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O.1002525E*05 U-5657 S68E*02 


value when a ray travels from 10 to 100 meters depth. The 
Significance of these three tables is that the ray acoustics 
methods are yielding the same phase values hy as the direct 
integration method. 

The conclusion to be derived from this is that ray 
acoustics theory can be used to accurately evaluate the phase 


integral. 


C. VARYING SOUND-SPEED GRADIENTS 
Both ray acoustics methods will now be applied to three 
adifferent test cases. A zero gradient, positive gradient 
(+0.016 sec), and a negative gradient (-0.016 sec’) linear 
sound-speed profile will be used. These simulations will show 
that both methods produce equivalent results within an 
acceptable tolerance for a variety of media. 
The conditions common to each simulation run are: 
- source depth of 10 meters, 
- horizontal range of 1 kilometer, 
- ocean depth of 100 meters, 


- speed of sound at the surface of 1500 meters per 
second, and 


- five values of radial spatial frequency (FR) 
evaluated. 
The results are presented in Tables 3.4 through 3.9. The 
phase integral solution is the phase change in the Y (depth) 


direction listed as the MODULO TWOPI THETAY value in radians. 
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TABLE 3.4 
ZERO GRADIENT, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 


TOTAL CPU TIME = 15 MIN . 52.38 SEC 


ey 


CASE 1A 
NDATA =» 4 NUM8ER OF GRADIENTS =» 3 
GC 11 ® O.00000EF*00 I/SEC Gt 2) = 0.00000F*00 I/SEC 
GC 31 = 0.00000€*00 I/SEC 
YO « 10.0 M DEPTH » 1300.0 M FREOC =» 250.0 HZ CMAX »® 3500.0 M/SEC 
FREQC/CMAX = 0.16667 CYCLES/M RATIO =» 0.9999 FRMAX »® 0.16665 CYCLES/M NFRS ® § 
DLTFR =» 0.033330 CYCLES/M RNGSTP »® 1.04 HRZRNG » 1000.0 M 
NOTE: FRMIN = DLTFR 
NFR BETAO(DEG) DEPTM(M) TRYLTC(SEC) SETAY(DEG] THETATCRAD] THETAR (RAD) THETAYC(RAD] MODULO TROP] 
THETAY CRAD) 
I 11.836 90.5} 3.333667 168.464 O.S236511E*04 0.2094186E°03 0.5027093E*046 0.544598} 
2 23.576 98.44 1.666833 1$6.426 0.2618256F°04 O.4188371E°03 0.21399619F°06 0.3037135 
3 36.866 56.46 Re121222 143.136 0.1745506E*04 0.6282557E*03 0.1117248E&*04 5.124306} 
4 $3.122 39.79 0.833417 126.878 0.1309128E*04 0.83767463E*03 0.4714836E203 0.2146886 
S 89.190 26.24 0.666733 89.190 O.10G7302E*04 0.10467093E*06 0.2094500E*00 0.2094500 
TOTAL CPU TIME &© O MIN . 22.63 SEC 
TABLE 3.5 
ZERO GRADIENT, ODE SOLVER SOLUTION 
RAY TRACING USING AKIMA CUBIC SPLINE & ODE SOLVER 
CASE 2A 
NDATA ® 4 
YO = 10.0 M DEPTH ® 100.0 mM FREQC = 250.0 HZ CMAX ® 1500.0 M/SEC 
FREQC/CMAX = 0.16667 CYCLES/M RAT1O & 0.9999 FRMAX »® 0.16665 CYCLES/M NFRS = 5 
DLTFR = 0.033330 CYCLES/M RNGSTP =» 1.04 HRZRNG 1000.0 M 
NOTE: FRMIN ® DLTFR 
NFR SETAOC( DEG) DEPTH{M) TRYLT(SEC) 8ETAY( DEG) THETATCRAD) THETAR (RAD) THETAY(RAD) MODULO THOPI 
THETAYCRAD) 
1 11.536 90.51 3.333667 168.464 0.S236S511E*04 0.20946186E*03 0.5027093E+04 0.544595} 
2 23.876 98.46 1.666833 156.424 0.2618256E°064 O.4188371E*03 0.2199419E 204 0.3037135 
3 36.866 56.46 B.211222 143.136 0.1745$04E*04 0.62825$7E°03 0.11317248E*04 5.124306} 
4 53.122 39.79 0.833617 126.978 O.2309128E*04 0.83767463E*03 0.4714536E°038 0.2146886 
S 89.190 26.16 0.666733 89.190 0.1047302E*04 0.10467093E*04 0.2094500F°00 0.2094500 


TABLE 3 


oO 


POSITIVE GRADIENT, PIECEWISE LINEAR SOLUTION 
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RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 18 
NDATA » 4 NUMBER OF GRADIENTS «© $3 
Gt 1) =# 0O.16000E-O1 I/SEC Gt 2) 2 0.16000E-01 I/SEC 
Gt 3) » 0.16000E-01 I/SEC 
YO : 10.0 DEPTH =» 100.0 FREOC # 250.0 H2 CMAX =» 1501.6 M/SEC 
FREOC/CMAX #® 0,1666° CYCLES/M RATIO = 0.9999 FRMAX » 0.16647 CYCLES/M NFRS ®&® § 
DLTFR » 0.033296 CYCLES/M RNGSTIP 2© 41.048 MRERNG «= 1000.0 
NOTE: FRMIN © DLTFR 
NFR BETAO(DEG) DEP TH(M) TRVLT(SEC) SETAY( DEG) THETAT(RAD) THE TAR (RAD) TMETAY(RAD) MODULO THOP] 
THE TAY(RAD) 
1 11.825 87.80 3.333652 168.466 O.S23648B8E*06 0.20919546E*03 0.5027293E +04 0.74662569 
2 23.552 96.99 1.666826 156.4625 O.26182¢61E*046 0.¢46183908E*03 0.21°98S50E +04 0.73468999 
3 36. 826 55.38 V.d221¢2 143.155 0.1765656E*04 0.6275863E*03 O.1117870E°0¢ §.7462108 
"6 53.0469 38.73 0.833378 126.927 O.1309066EF*04 O.8367817E*03 0.4722847E*03 1.0657792 
5 87.3683 50.71 0.667000 87.974 O.1067721E*06 0.1065977E*06 0.1743798E*01 1.7437978 
TOTAL CPU TIME » O MIN. 22.05 SEC 
TABLE 3.7 
POSITIVE GRADIENT, ODE SOLVER SOLUTION 
RAY TRACING USING AKIMA CUBIC SPLINE 8&8 ODE SOLVER 
CASE 28 
NDATA = 4 
YO » 10.0H DEP TM »& 100.0 FREOC = 250.0 #2 CMAX = 1501.6 M/SEC 
FREQC/CMAX # 0.1664% CYCLES/M RAT1O =» 0.9999 FRHAX = 0.16647 CYCLES/H NFRS = 5§ 
OLTFR »&= 0.033294 CYCLES/H RNGSTP » 1.04 HRZRNG » 1000.0 4 
NOTE: FRMIN ©» DLTFR 
NFR BETAO( DEG) DEPTH(M) TRVLT(SEC) SETAY( DEG) THETAT(RAD) THE TAR (RAD) THETAY(RAD) MODULO THOP] 
THE TAYC(RAD) 
1 11.525 87.80 $.333652 168.466 0.5236488E°04 0.2091956E*03 0.5027293E*04 0.7442693 
2 23.552 96.99 1.668826 156.425 O.26182461EFE*04 O.4183908E*03 0.21°°8S0E*04 0.7349005 
3 36.8624 55.38 L.2221292 363.155 0.17465456E*04 0.627586S5E*03 O.1117870E+04 §.74662125 
4 53.049 "8.73 0.852378 126.927 O.1309066E*0% O.8367817E*03 0.46722847E*03 1.0657783 
Ss 87.363 §0.71 0.667000 87.974 0.1067721E*04 0.1065977E*04 0.1743798E*0)1 1.7437978 
TOTAL CPU TIME © 16 MIN . 24.30 SEC 


=? 


TABLE 3.8 


NEGATIVE GRADIENT, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPRDXIMATION 
CASE 1C 


NDATA = 4G NUMBER OF GRADIENTS ® 3 
GO 1) = -0.16000E-01 1/SEC GO 2) = -0.16000E-01 1/SEC 
GC 3) 2 -0.16000E-01 1/SEC 
Yo = 10.0 M DEPTM = 100.0 M FREOC =» 250.0 HZ CMAX = 1500.0 M/SEC 
FREQC/CMAX = 0.16667 CYCLES/M RATIO = 0.9999 FRMAX = 0.16665 CYCLES/M NFRS = § 
DLTFR = 0.033330 CYCLES/M RNGSTP = 1.0 HRZRNG # 1000.0 M 
NOTE: FRMIN #® DLTFR 
NFER BETAO(DEG) DEPTM(M) TRYLT(SEC) SETAY(DEG) THETAT(CRAD) THETAR(RAD) THETAYC(RAD) MODULD TWOP) 
THETAY (RAD) 
1 11.535 87.78 2.337240 168.475 0.5242126E*04 0.2094186E*03 0.50352706E°04 6.1576516 
2 25.575 96.98 1.668622 156.450 0.2621065E*04 0.46188S71E*0$ 0.2202228E°04 3.1134653 
3 36.861 $5.32 1.112478 143.160 0.17676146E204 9 0.6282557E*0S 0.1119159E°04 G. 7515715 
“ §3.114 38.63 0.834245 126.909 0.1310586E*04 0.8376743E*03 0.4729117E*08 1.672792) 
+) 68.835 35.67 0.667046 88.224 0.1047793E*04 0.1047093E*04 0.6999612E*00 0.6999612 
TDTAL CPU TIME = 0 MIN » 22.61 SEC 
TABLE 3.9 


NEGATIVE GRADIENT, ODE SOLVER SOLUTION 


RAY TRACING USING AKIMA CUBIC SPLINE & ODE SOLVER 
CASE 2C 


FREQC = 250.0 HZ 


CMAX = 


1500.0 M/SEC 


NDATA 2 4 

yo «= 10.0 M DEPTH = 100.0 M 

FREQC/CMAX = 0.16667 CYCLES/M RATID = 0.9999 

DLTFR = 0.033330 CYCLES/M RNGSTP = 1.06 HRZRNG « 

NOTE: FRMIN = DLTFR 

NFR BETAO(DEG) DEPTH(M) TRYLTC(SEC) SETAY(DEG) THETATCRAD) 
1 11.535 87.77 3.337240 168.475 0.52421246E +04 
2 23.578 96.98 1.668622 156.450 0.2621065E+04 
3 36.861 55.32 1.112688 143.160 0.1747414E*04 
4 S$3.114 38.63 0.834345 126.909 0.1310586E+*04 
S 88.835 35.67 0.667046 88.224 0.1047792E*04 


TOTAL CPU TIME = 16 


MIN » 21.26 SEC 


ee, 


1000.0 M 


THE TAR CRAD} 


0.2094186E*03 
0.4188S71E°03 
0.6282557E*03 
0.8376743E*03 
0.104709S5E*04 


FRMAX = 0.16665 CYCLES/M NFRS = 5 


THETAY(RAD} 


0.5032706E*04 
0.22020228E*04 
0.11191S9E*04 
0.46729117E*03 
0.6999612E*00 


MDDULD TWOP 1 

THE TAYC(RAD) 
6.1576656 
3$.1134683 
0.7515729 
1.6727921 
0.6999612 


Tables 3.4 and 3.5 show the results for the zero gcadvenas 
case. A comparison between two tables shows that the phase 
values (MODULO TWOPI THETAY) agree perfectly. Additionally, 
all other calculated values in the tables agree. The only 
difference is the TOTAL CPU TIME. Both methods are using a 
conservative range step (RNGSTP) of 1 meter, i.e., the 
incremental increase in range is 1 meter in the propagation 
calculations. For these conditions, the ODE solver is taking 
over 30 times longer to run the simulation. 

Tables 3.6 and 3.7 show the results for the positive 
gradient case. A comparison between these two tables, as 
well as the negative gradient results Shown in Tables 3.8 and 
3.9, shows acceptable agreement between the two very 
different calculation methods. The only discrepancies in 
calculated values occurs in the modulo 27 0, phase values 
needed for solving the phase integral. The phase calculation 
is most challenging because of the nature of phase (being a 
modulo 27 function). The discrepancies for both the positive 
and negative gradients are on the order of hundred-thousandths 
of a radian or less. The CPU times continue to follow the 
pattern seen earlier that the ODE solver requires over 30 
times longer to complete a simulation run. 

The following ray trace plots, Figures 3.1 through 3.6, 
correspond to Tables 3.4 through 3.9. The plots assist 


interpretation of the data, but as expected from the tabular 
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RAY TRACE: AKIMA CUBIC SPLINE & ODE SOLVER 
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Figure 3.5 
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Ray trace corresponding to Table 3.9 


Figure 3.6 
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results, no variations in depth or angle of arrival can be 


detected. 


D. CHALLENGING TEST CASES 

A smooth sound-speed profile in the form of a half 
Sinusoid describes the medium for the next set of 
Simulations. To this point, the agreement in accuracy between 
the two methods has been excellent; however, the shape of a 
Sine curve is difficult to approximate with piecewise linear 
segments. While this will be a challenging test, smoothly 
curved profiles do occur naturally. The following set of 
values were used: 


- 1500 m/sec is the sound speed at the ocean surface and 
at the 2000 meter bottom, and 


- a minimum sound speed of 1475 m/sec occurs at the 1000 
meter depth. 


Three sets of simulations are presented using 5, 11, and 
17 sound-speed versus depth data pairs. The simulation 
results will be examined to quantify the number of ocean 
medium sound-speed samples required by each method to converge 
to a solution. Also, the agreement in results for the two 
methods and the CPU times will be evaluated. 
As the number of equally spaced data pairs varies, the 
medium and source conditions will be: 
- source depth of 1000 meters, 
- horizontal range of 3 kilometers, and 
- five values of radial spatial frequency (FR) 


evaluated. 
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Tables 3.10, 3.11, and 3.42 show the results fon thes 
11, and 17 data pair test cases for the piecewise linear 
method. In comparing the first two of these tables, large 
discrepancies are seen in the final depths for each of the 
five radial frequencies evaluated—errors as large as seven 
meters. Likewise, travel time values differ by hundredths of 
a second meaning that modulo 2m” phase values’ show no 
correlation between test cases. Unfortunately, the variations 
between Tables 3.11 and 3.12 where more data points are used, 
show the piecewise linear method is not converging to a 
Ssoluteron. 

Figures 3.7 through 3.9 correspond to Tables 3.10 through 
3.12. While numerical phase values show gross disagreement, 
the inability to detect differences in the graphical 
representations illustrates the very sensitive nature of the 
phase calculations. 

Tables 3.13 through 3.15 show the results when the ODE 
solver runs the simulations for 5, 11, and 17 data pairs. The 
5 data pair case of Table 3.13 again shows significant 
disagreement with the Table 3.14 results using 11 data pairs. 
Final depths show differences of as much as a meter while 
travel times vary by milliseconds. As with the piecewise 
linear method, the modulo 2m” phase values have not yet 
converged to a solution. Variations of up to 2 radians are 


seen. 
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TABLE 3.10 


FIVE DATA PAIRS, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 


CASE 13805 
NDATA = § NUMBER OF GRADIENTS = 4 
GC 1) = ~0.35355E-01 1/SEC Gt 2) = ~0.14646S5E-01 1/SEC 
GC 3) 3 0.146465E-01 1/SEC G( 4) = 9Q.3S3SSE-01 1/SEC 
ERAYS 3 F PRTALL = F 
YO = 1000.0 M DEPTH © 2000.0 M FREQC = 250.0 HZ CMAX = 1500.0 M/SEC 
FREQC/CMAX = 0.16667 CYCLES/M RATIO = 0.9999 FRMAX = 0.16665 CYCLES/M NFRS = § 
DLTFR = 0.033330 CYCLES/M RNGSTP = 1.0 M HRZRNG 2 3000.0 M YR = 40.0 M 
NOTE: FRMIN 2 DLTFR 
NFR BETAOQ(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETATCRAD) THE TAR (RAD) THETAY(RAD) MODULO TWOP] 
THE TAY(RAD) 
1 11.341 144.74 10.207186 168.504 0.1603341E*05 0.6282557E*03 0.1540515E*05 5 .0674943 
2 23.160 41.81 5.103311 156.449 0.80162626*04 90.1256511E*04 0.6759751E°04 § 32646822 
3 36.183 1064.46 3.402439 36.180 0.834646S38E*04 0.1884767E*04 0.3459771E*04 4°.0191220 
4 51.868 682.47 2.854747 127.902 0.4012987E*04 0.2513023E°04 0.1499964E*04 4.5658515 
5 79.494 1510.35 2.057830 81.249 O0.32324633E*04 0.3141278E*04 0.9115404E*02 3.18 94409 
TOTAL CPU TIME = 0 MIN » 24.95 SEC 


ELEVEN DATA PAIRS, 


PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 


CASE 181) 

NDATA = 11 NUMBER OF GRADIENTS * 10 

GC 1) 2 ~0.38627E-01 I1/SEC G( 2) = -0.348466E-01 1/SEC 

GC 3) 2 -0.27654E-01 1/SEC GC 4) s -0.17755E-01 1/SEC 

GC 5) = -0.61179E-02 1/SEC Gt 6) = 0.61179E-02 1/SEC 

G( 7) = 0.1775S5E-01 1/SEC G( 8) = 0.27654E-01 1/SEC 

GC 9) = 0.34846E-01 1/SEC G(10) = 0.38627E-01 1/SEC 
ERAYS = F PRTALL = F 

YO = 1000.0 M DEPTH = 2000.0 M FREOC = 


FREOC/CMAX = 0.16667 CYCLES/M 
0.033330 CYCLES/M 


DLTFR = 
NOTE: FRMIN = DLYFR 


NFR BETAO(DEG) DEPTH(M) TRYLT(SEC) BETAY(DEG) 


1 11.341 137.50 
r4 23.160 37.97 
3 36.153 1067.39 
4 51.868 679.33 
5 79.494 1§20.40 


10.216764 
5.108073 
3.405626 
2.557346 
2.060168 


TOTAL CPU TIME » 0 MIN , 25.36 SEC 


RATIO = 0.9999 
RNGSTP 2 


1.0 M 


THETATC(RAD) 


168.506 0.1604846E°05 
1$6.449 0.8023742E +04 
36.165 0.5349545E*04 
127.965 0.4017070E+04 
81.354 0.32361046E +04 
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HRZRNG 2 


0.6282557E*03 
O.1°056S11E*04 
0.1884767E*04 
0.2513023E*04 
0.3141278E*04 


250.0 H2 CMAX 8 
FRMAX = 0.16665 CYCLES/M NFRS = § 
3000.0 M YR = 40. 
THETAR (RAD) THE TAY (RAD) 


0.1542020E*05 
0.6767231E*04 
0.3464778E +04 
0.1504068E+04 
0.9482580E*02 


1500.0 M/SEC 


0 M 


MODULO TWOP 1 
THE TAY( RAD) 
1.2628446 
0.2400057 
2.7430820 
2.366369) 
0.5780165 


TABLE 3a 12 
SEVENTEEN DATA PAIRS, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 


CASE 13817 

NDATA 8 17 NUMBER OF GRADIENTS = 16 

GC 1) = -0.39018E-01 1/SEC GC 2) # -0.37519E-01 17SEC 

G( 3) 2 ~0.34577E-01 1/SEC GC 4) © -0.30307E-01 1/SEC 

G( 5) = -0.24873E-01 1/SEC GC 6) = ~0.18482E-0] 1/SEC 

G( 7) = -0.11381E-01 1/SEC G( 8) = -0.38429E-02 1/SEC 

GC 9) = 0.38429E-02 1/SEC G(10) = O.11381E-01 1/SEC 

GO1ll) = 0.18482E-01 1/SEC G12) = 0.24873E-01 I/SEC 

G(13) = 0.30307E-01 1/SEC G16) = 0.34577E-01 1/SEC 

G(15) = 0.37519E-01 1/SEC G(16) = 0.39018E-01 1/SEC 

ERAYS 3 F PRTALL = F 

YO = 1000.0 M DEPTH = 2000.0 M FREOQC s 250.0 HZ CMAX = 1500.0 M/SEC 

FREQC/CMAX = 0.16667 CYCLES/M RATIO = 0.9999 FRMAX = 0.16665 CYCLES/M NFRS #» 5 

DLTFR = 0.033330 CYCLES/M RNGSTP = 1.04 HRZRNG 2 3000.0 M YR = 40.0 M 

NOTE: FRMIN = DLTFR 

NFR BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETAT(RAD) THE TAR CRAD ) THETAY(RAD) MODULO TWOP! 

THETAY (RAD) 

| 11.341 136.67 10.217866 168.506 O.160S019E*0S 0.6282557E*03 0.1542193E+*05 2.9937577 
2 23.160 37.53 §.108620 156.449 0.8024602E°04 O.1256511E*04 0.6768091E*04 1.1002238 
3 36.153 1067.73 3.405996 36.160 0.5350126E*04 0.1884767E*04 0.3465359E+04 3.3242727 
4 51.868 678.98 2.557644 127.973 0.4017538E*04 0.2513023E*04 0.1504515E*04 2.8338755 
> 79.494 1521.62 2.060452 81.322 0.32365S0E*04 0.3141278E*04 0.9527150E*02 1.0237247 


TOTAL CPU TIME = 0 MIN » 25.94 SEC 
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FIVE DATA PAIRS, ODE SOLVER 


TABLE 3.13 


RAY TRACING USING AKIMA CU8SIC SPLINE & ODE SOLVER 


CASE 2505 


NDATA = § 


ERAYS = F 
YO = 1000.0 M 


FREOC/CMAX = 0.16667 CYCLES/H 


DLTFR » 
NOTE: FRHIN = 


NFR BETAQ(DEG) DEPTM(M) TRYLT(SEC) BETAY(DEGI 


1 11.34] 
2 23.160 
3 36.153 
4 51.868 
5 79.494 


DLTFR 


134.85 
36.56 
1068.49 
678.30 
1521.82 


0.033330 CYCLES/M 


PRTALL © F 


DEPTH = 2000.0 M 
RATIO = 0.9999 


FREOC = 


250.0 HZ CMAX = 15 

FRMAX = 0.16665 CYCLES/H NFRS = § 
3000.0 M YR «= <0. 

THE TAR (RADI THETAY (RADI 
0.62825S7E*03 0.1842569E*05 
O.12565116°06 0.6769980E*04 
0.18847676°04 0.3466635E°04 
0.25130236°06 0.1505379E*04 
0.31412786°04 0.95344668E*02 


RNGSTP = 1.046 HRZRNG = 
THETAT(RAD I 

10,220287 168.509 0.1605394E*05 
§.10°823 156.452 0.80264691E6°046 
3.406608 36.157 0.535146006°04 
2.558194 127.973 0.4018402E°04 
2.060498 81.306 0.32366236°04 

TABLE 3.14 


ELEVEN DATA PAIRS, ODE SOLVER 


RAY TRACING USING AKIMA CUBIC SPLINE & ODE SOLVER 


00.0 M/SEC 


0M 


MODULO THOPI 
THETAY (RAD) 
- 0.4666177 
2.98918646 
4.59970158 
$.6978733 
1.0968988 


CASE 2$11 
NDATA = 1] 
ERAYS » F PRTALL © F 
YO = 1000,0 M DEPTH = 2000.0 M FREOC = 250.0 H2 CMAX = 1500.0 M/SEC 
FREOC/CHAR = 0.16667 CYCLES/M RATIO = 0.9999 FRMAX = 0.16665 CYCLES/M WFRS = 5 
OLTFR © 0.033330 CYCLES/M RNGSTP = 1.04 HRZRNG ® 3000.0 M YR = 40.0 M 
NOTE: FRMIN © DLTFR 
NFR BETAO(DEG!] DEPTHI(M!] TRYLT(SEC) BETAY(DEG] THETAT(RAD) THE TAR (RAD) THETAY(RAD) MODULO THOP] 
THETAY(RAD) 
} 11.341 136.09 10.216623 168.806 O.1605138E°0S$ 0.62825S7E°03 0.15462312E°05 4.1828982 
2 23.160 37.23 §$.108997 156.449 0.8025194E°04 O.1256511E°06 0.6768685E°04 1.6920625 
3 36.153 1067.97 3.406254 36.157 O.S3S0SS1E°04 0.16846767E°04 0.34687646E*04 3.728877) 
4 §1.868 678.73 2.$87668 127.978 0.40178S8E*06 0.2513023E°04 0.150483SE*04 3.1540595 
§ 79.496 1822.39 2.060633 81.312 O.S2368346E204 0.314612786°04 0.985S5SS59E902 1.3078080 
TOTAL CPU TIHE = 48 HIN . 16.83 SEC 
TABLE 3.15 


SEVENTEEN DATA PAIRS, 


RAY TRACING USING AKIHA CUBIC SPLINE & ODE SOLVER 
CASE 2317 


NDATA = 17 


ERAYS = F 


Yo 


FREOC/CHAX = 0.16667 CYCLES/M 
DLTFR = 0.033330 CYCLES/H 


NOTE: FRMIN = 


NFR BETAOCDEGI 


® 1000.0 M 


DL 


TFR 


1 11.341 136.13 
2 23.160 37.25 
3 36.183 1067.96 
4 51.868 678.74 
5 79.494 1822.36 


PRTA 


DEPTM = 2000.0 M 
RATIO = 0.9999 
RNGSTP s 


10.216580 1 
$.108975 1 
3.406239 
2.557636 1 
2.060625 


TOTAL CPU TIME = 47 MIN , 45.20 SEC 


LL © F 


DEPTH(MI TRVLT(SEC) BETAY(DEG) 


68.506 
56.469 
36.157 
27.978 
81.317 


ODE SOLVER 


FREQC = 250.0 M2 CMAX = 

FRHAX © 0.16665 CYCLES/M NFRS = § 
La] KRZRNG = 3000.0 M YR = 
THETAT(CRAD) THETAR (RADI THE TAY TRAD) 


0.160S131E°05 
0.8025159E*04 
0.§350S07E*04 
0.4017839E°04 
0.3236822E°04 
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0.62825S57E°03 
0.12565 116 2046 
0.18684767E°04 
0.25130236°04 
0.3141278E*04 


0.15462305E°05 
0.6768648E *04 
0.24657406 +046 
0.1504817E*0<¢ 
0.9554374E*02 


1§00.0 M/SEC 


40.0 MK 


MODULO TROP] 
THETAY(RAD 1 


4.1156876 
1.6570617 
3.7048568 
3.1353514 
1.2959589 


A comparison of Tables 3.14 and 3.15 shows that phase 
calculation results converged for each of the five FR values. 
The increase in data pairs to 17 from 11 has changed the 
final depth values by only a few hundredths of a meter. 
Likewise, the very sensitive modulo 2T phase values show a 
change of only hundredths of a radian or less. The ray trace 
plots for Tables 3.13 through 3.15 are provided in Figures 
feeO through 3.12. 

The CPU times for these three tables show the ODE solver 
to be costly, but insensitive to increases in the number of 
data pairs. That 1S, approximately the same amount of CPU 
time 1S required regardless of the number of data pairs used. 
The CPU time for the piecewise linear method also proved 
fairly insensitive to the number of data pairs used as seen in 
Tables 3.10 through 3.12. 

The significant findings of this section are: 

- the ODE solver can perform accurate phase 

calculations with only 17 data pairs sampled from a 

2000 meter deep ocean having a curved, sinusoidal 


sound-speed profile, 


- the piecewise linear approach does not converge to a 
solution using the 17 data pairs, and 


- the differential equation solver is very costly to use 
in terms of CPU time versus the piecewise linear approach. 
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Figure 3.10 
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Figure 3.11 
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E. SPECIAL INVESTIGATIONS 

Each of the two phase calculation methods has shown one 
characteristic weakness. The piecewise linear method did not 
converge to a solution in the sinusoidal sound-speed profile 
case when 17 data pairs or samples were provided. In 
contrast, the ODE solver does converge to a solution, but is 
costly in terms of CPU time. This section will attempt to 
tailor the settings for each method to compensate for its 
weaknesses. 

The simulation results in Tables 3.16 through 3.18 were 
produced by the piecewise linear method using 29, 55, and 65 
data pairs, respectively. All medium and sound-speed profile 
characteristics remain the same as in the previous section. 
As the number of data pairs increases, the 0, values are 
converging to the Table 3.15 solutions obtained with the ODE 
solver. Unfortunately, the phase values continue to show 
unacceptable errors of tenths and hundredths of a radian for 
the best case of 65 data pairs. 

Table 3.19 shows the results for the ODE solver when the 
range step (RNGSTP) parameter is increased to five meters. As 
explained earlier, the range step is the differential 
equation system independent variable. All calculations prior 
to Table 3.19 used a conservative range step of one meter for 
both methods. Table 3.19 compares favorably with Table 3.15 
achieving a balance between CPU time and accuracy. While the 


degradation in accuracy is only microseconds for travel time 
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TABLE 3.16 
TWENTY-NINE DATA PAIRS, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1829 


NDATA #® 29 NUMBER OF GRADIENTS ® 28 

GC 1) # -0.39188E-01 I/SEC G( 2) = -0.38695E-01 1/SEC 
G( 3) # -0.37715E-01 I/SEC G( 4) = -0.36262E-01 1/SEC 
G( 5) # -0.34352E-01 I1/SEC GC 6) # -0.32010E-01 1/SEC 
GC 7) # -0.29266E-01 I/SEC G( 8) = -0.26154E-01 1/SEC 
GC 9) 8 -~0.22712E-01 1/SEC G(10) = -0.18986E-01 1/SEC 
G(11) © -0.15020E-01 I/SEC G(12) # -0.10866E-01 1/SEC 
G(13) © -0.65745E-02 1/SEC G(14) » -0.22007E-02 1/SEC 
G(15) #8 0.22007E-02 1/SEC G(16) = 0.65745E-02 1/SEC 
G(17) © 0.10866E-01 1/SEC G(18) © 0.15020E-01 1/SEC 
G(19) = 0.18986E-01 1/SEC G(20) # 0.22712E-01 1/SEC 
G(21) = 0.26154E-01 1/SEC G(22) # 0.29266E-01 I/SEC 
G(23) © O0.32010E-01 1/SEC G(24) 8 0.34352E-01 1/SEC 
G6(25) = 0.36262E-01 1/SEC G(26) = 0.37715E-01 1/SEC 
G(27) = G.38695E-01 1/SEC G(28) = 0.39188E-01 1/SEC 


ERAYS © F PRTALL = F 

YO = 1000.0 M DEPTH = 2000.0 M FREQC es 250.0 HZ CMHAX s 1500.0 M/SEC 
FREQC/CMAX = 0.16667 CYCLES/M RATIO = 0.9999 FRMAX = 0.16665 CYCLES/M WFRS = § 

DLTFR = 0.033330 CYCLES/H RNGSTP = 1.046 HRZRNG ® 3000.0 YR = 40.0 M 


NOTE: FRMIN © DLTFR 


NFR BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETAT(RAD) THETAR (RAD) THETAY(RAD) MODULO TWOP] 


‘ 


THETAY (RAD) 
11.341 136.31 10.218342 168.506 O.1605093E°0S 0.6282557E°03 0.1542268E*05 3.7413092 
23.160 37.34 §.108857 156.449 0.8024973E204 O.1256511E204 0.67684562E404 1.4710844 


36.153 1067.88 3.406158 36.157 0.5350380E*04 0.1884767E°04 0.3465613E+04 3.5775983 
51.868 678.82 2.557772 127.976 0.4017739E*04 O.2513023E°904 0.1504716E404 8.0348743 
79.494 1§22.11 2.060567 81.320 0.3236730E204 0.31241278E*04 0.9545197E*02 1.2041934 


wn ww AH = 


TOTAL CPU TIME = 0 MIN » 26.46 SEC 
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TABEGE, 3.17 
FIFTY-FIVE DATA PAIRS, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 


CASE 1S55 

NDATA = 55 NUMBER OF GRADIENTS = 84 

GC 1) # -0.59248E-01 1/SEC GC 2) = -O.39115E-01 1/SEC 

GC 3) = -0.388S0E-01 1/SEC G( 4) =» -0.38453E-01 1/SEC 

G( 5) = ~0.37926E-01 1/SEC G( 6) 2 ~0.37271E-01 1/SEC 

GC 7) 8 ~0.36490E-01 1/SEC G( 8) = -0.35586E-01 1/SEC 

GC 9) 2 -0.34561E-01 1/SEC G(10) = -0.33419E-01 I/SEC 

GC11) = -0.32163E-01 1/SEC G(12) = -0.30800E-01 1/SEC 

GC1S) = -0.29331E-01 1/SEC G(14) = -0.27764E-01 1/SEC 

G(1S) = -0.26103E-01 1/SEC G(16) = -0.2435535E-01 1/SEC 

G17) = -0.22521E-01 1/SEC G(18) = -0.20613E-01 1/SEC 

G(19) = ~-0.1865SE-01 1/SEC G(20) = -0.16594E-01 1/SEC 

G(21) = -0.14497E-01 1/SEC G(22) = -0.12350E-01 1/SEC 

G(23) = -0.10162E-01 1/SEC G(24) = -0.79399E-02 1/SEC 

G(25) = ~0.56907E-02 1/SEC G(26) #8 -0.34221E-02 1/SEC 

G(27) 8 -0.11420E-02 1/SEC G(28) &© O0O.11420E-02 1/SEC 

G(29) = 0.34221E-02 1/SEC G(30) = 0.56907E-02 1/SEC 

G(3$1) 8 0.79399E-02 1/SEC G(3$2) 8 0.10162E-01 1/SEC 

G(3$3) = 0.12350E~-01 1/SEC G(34) & 0.14497E-01 1/SEC 

G(35) = 0.16594E-01 1/SEC G(36) = 0.1863S5E-01 1/SEC 

G(37) = 6.20613E-01 1/SEC G(38) 2 0.22521E-01 1/SEC 

G(39) s 0.24353E-01 1/SEC G(40) = 0.26103E-01 1/SEC 

G(41) = 0.27764E-01 1/SEC G(42) 8 0O.29331E-01 I/SEC 

G(43) © O.30800E-01 1/SEC G44) © 0.32163E-01 1/SEC 

G(43) © 0.33419E-01 1/SEC G46) © 0.346561E-01 1/SEC 

G(47) = 0.35586E-01 1/SEC G(48) © 0.36490E-01 1/SEC 

G(49) = 0.37271E-01 1/SEC G(S0) = 0.37926E-01 1/SEC 

G(S1) = 0.38453E-01 I/SEC G52) = 0.38850E-01 1/SEC 

G(53) = O.3S9115E-01 I/SEC G(S4) © 0.39248E-01 I1/SEC 

ERAYS = F PRTALL = F 

YO = 1000.0 # DEPTH = 2000.0 M FREQC = 250.0 HZ CMAX = 1500.6 M/SEC 

FREQC/CMAX # 0.16667 CYCLES/M RAT1O = 0.9999 FRMAX #® 0.16665 CYCLES/M NFRS = § 

DLTFR = 0.033330 CYCLES/M RNGSTP # 1.0 M HRZRNG 8 3000.0 M YR 2 40.0 M 

NOTE: FRMIN ® DUTFR 

NFR BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETAT(RAD) THETAR (RAD ) THETAY(RAD) MODULO TNOPI 

THETAY (RAD ) 

1 11.341 136.18 10.218510 168.506 0.1605120E705 0.62825S57E503 0.1542294E08 4.0060532 
2 23.160 37.28 5.108940 156.449 0.8025104E*04 O.1256511E404 0.6768593E404 1.60247883 
3 36.153 1067.94 3.406215 3¢.157 0.8350470E°04 0.1884767E404 0.3465702E404. 3.6673828 
& 51.868 678.77 2.857817 127.977 0.4017810E*04 0.2513023E*04 0.1504787E+04 3.1058334 
5 79.494 1522.29 2.060608 81.317 0.3236795E*04 0.3141278E704 0.9551684E%02 1.2690645 


TOTAL CPU TIME # 0 MIN » 26.61 SEC 
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TABLE 3.18 
SIXTY-FIVE DATA PAIRS, PIECEWISE LINEAR SOLUTION 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATION 
CASE 1565 


NDATA = 65 NUMBER OF GRADIENTS #= 64 

GC 1) = -0.39254E-0] I/SEC GC 2) = -0.39160E-01 1/SEC 
GC 3) = -0.38971E-01 1/SEC GC 4) = -0.38688E-01 1/5EC 
GC 5) s -0.38312E-01 I1/SEC GC 6) = -0.3784G6E-01 I1/SEC 
G( 7) = -0.37284E-01 1/SEC G( 8) = -0.56635E-01 I/SEC 
GC 9) = -0.35897E-01 1/SEC G(10) « -0.35073E-01 1/SEC 
GC11) = -0.34165E-01 1/SEC GCl2) = -0.33174E-01 I/SEC 
G(13) = -O.32103E-01 I/SEC G(14) = -0.30955E-01 I/SEC 
GO1S) © -0.29733E-01 1/SEC G(16) = -0.28438E-01 I/SEC 
G(17) s -0.27075E-01 I1/S5EC G(18) 5 -0.25647E-01 1/SEC 
G(19) =» -0.24158E-01 I1/SEC G(20) = -0.2¢2610E-01 I1/SEC 
G21) = -0.21007E-01 I/SEC Glco) * -0.1935GE-01 1/SEC 
G(235) = -0.17654E-01 I/SEC G(24) = -0.15912E-01 I/SEC 
G25) # -0.j}61352E-0) 1/SEC Gll6) £ -OU 12S 17E-C1 17SEC 
G(27) = -0.10473E-01 I/SEC G(28) = -0.86032E-02 I/SEC 
G(29) = -0.67130E-62 I/SEC G(30) = -0.48066E-02 I/SEC 
G(31) = -0.28886E-02 1/SEC G(32) = -0.9636G4E-03 I/SEC 
G(33) = 0.96364E-03 I/SEC G(34) = 0.28886E-02 1/SEC 
G(35) = 0.48066E-02 I/SEC G(36) = O.67130E-02 I/SEC 
G(37) = 0.86032E-02 1/SEC G(l8) = 0.10473E-01 1/SEC 
G(39) = 6.12317E-01 I1/SEC G(40) = O.14132E-601 I/SEC 
G(41) = 6.15912E-01 1/SEC G(42) 2 0.17654E-01 I1/SEC 
G(63) = 0.19354E-01 I/SEC G(44) 8 0.21007E-01 I/SEC 
G(465) = 0.22610E-01 I/SEC G(46) = 0.24158E-01 I/SEC 


G(47) = 0.25647E-0) I/SEC G(48) = 0.27075E-01 I/SEC 

G(49) = 0.28438E-01 I1/SEC G(50) = 0.29733E-01 I/SEC 

G(S1) = 6.309SSE-01 I1/SEC G(52) = O.32103E-01 I/SEC 

G(S3) = 6.33174E-0) 1/SEC G(54) = 0.34165E-01 I/SEC 

G(SS) = 0.35073E-01 1/SEC G(S6) = 06.38897E-01 1/5EC 

G(S57) = 0.26635E-01 1/SEC G(58) = 0.37284E-01 I1/SEC 

G(59) = 6.3784GE-01 1/S5EC G(60) = O.38312E-01 I/SEC 

G(61) = 0.38688E-01 1/SEC G(62) 2 0.38971E-0) I/SEC 

G(63) = 6.39160E-01 1/SEC G(64¢) = 0.39254E-0) I/SEC 

ERAYS = F PRYTALL = F 

YO » 1000.0 M DEPTH = 2000.0 M FREQC = 250.0 HZ CMAX = 1500.0 M/SEC 
FREQC/CMAX = 0.16667 CYCLES/M RATIO = 0.9999 FRMAX = 0.16665 CYCLES/M NFRS » 5 
DLTFR = 0.032330 CYCLES/M RNGSTP = 1.04 HRZRNG 2 3000.0 M YR = 40.0 


NOTE: FRMIN » DLTFR 


NFR BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETAT(RAD) THE TAR (RAD) THETAY(RAD) MODULO TWOP I 


THE TAY (RAD) 
11.34) 136.17 10.218528 168.506 0.1605123E*05 0.6282557E*0S 0.1542297E*05 4.0340889 
23.160 37.27 5.108949 156.449 O.80°05118E*04 O.1256511E*06 0.6768607E+04 1.6163983 


36.153 1067.94 3.406221 BIR evry 0.5350479E*0G 0.1884767E*0G 0.3465712E*04 3.6769515 
51.868 678.76 Z2.9978ca. Vel 978 0.4017818E*0G4 O.25123023E*04 0.1504795E*04 3.113625] 
79.494 1522.31 2.060612 81.317 0.3226802E*04 0.3141278E*04 0.9552360E*02 1.2758237 


1v ee > os 


TOTAL CPU TIME = 0 MIN », 26.84 SEC 
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and hundreds or thousandths of a radian for modulo 27 6,, the 
reduction in CPU time is following a nearly linear relation 
to the increase in range step, that is, CPU time was reduced 


approximately by one-fifth. 


F. 


ray acoustics methods. 
these two tables. 
be that 


because 


v 


TABLE 
INCREASED RANGE STEP, 


RAY TRACING USING AKIMA CUBIC SPLINE & ODE SOLVER 


CASE 2517 


NDATA = 17 


ERAYS = F 


YO = 1000.0 M 
FREQC/CMAX = 0.16667 CYCLES/M 
0.033330 CYCLES/M 


OLTFR = 


NOTE: FRMIN = DLTFR 


NFR BETAO(DEG) DEPTH(M) TRVYLT(SEC) BETAY(DEG) 


11.34} 
23.160 
36.153 
51.868 
79.496 


fe >) ee ee 


TOTAL CPU TIME 


136.09 
37225 
1067.96 
678.74 
1$22.36 


s 10 MIN. 3.53 SEC 


EIGENRAYS 


10. 
Ss. 
3. 
2. 
2. 


occurs in nature. 


meters depth and a local minima occurs at the ocean surface. 


The tables show that each method found the same eigenrays at 


PRTALL = F 


DEPTH = 2000.0 M 
RATIO = 0.9999 


RNGSTP = 


218584 
108975 
406239 
557836 
060625 


The figures show the sound-speed profile to 


The minimum sound-speed occurs at 1000 


168. 
156. 
36. 
127. 
81. 


$06 
449 
137 
978 
31? 


Figures 3.13 and’ 3.14 correspond to 


of a SOFAR channel. 


5.0 


o oa Oo 90 © 


HRZRNG 2 


THE TAT (RAD) 


©1605131E+05 
-8025159E +04 
-S3S0S07E +04 
-4017839E +04 


Sau, 
ODE SOLVER SOLUTION 


FREQC = 


250.0 HZ CMAX = 
FRMAX = 0.16665 CYCLES/M NFRS =» 5 
3000.0 ™ YR = 
THE TAR (RAD) THETAY (RAD) 
0.628255S7E*03 0.1542306E*05 
O.1256511E*04 0.6768648E+04 
0.1886767E*04 0.34657460E°0¢6 
0.2513023E*04 0.1504817E+04 
0.3141278E*04 0.9554379E*02 


-S236822E*06 


40.0 


wo = — kd kot 


1500.0 M/SEC 


MODULO TWOP 1 
THE TAY (RAD) 
4. 
» 6570457 
- 7047770 
- 1353500 
-2960153 


1219846 


Tables 3.20 and 3.21 show the eigenrays found by the two 


w+ 


This channel profile was chosen 
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it produces interesting ray patterns, and it commonly 


TABLE 3.20 
EIGENRAYS FOUND USING THE 
PIECEWISE LINEAR APPROXIMATION METHOD 


RAY TRACING USING PIECEWISE LIHEAR APPROXIMATION 


CASE I1CM4 
NDATA = 4 NUHBER OF GRADIENTS = 3 
Gt 1) = 0O.16000E-01 I/SEC GC 2) = -0.1888°E-01 I1/S5EC 
GC 3) = O.17000E-01 1/S5EC 
ERAYS = T PRTALL = F 
YO = 100.0 H DEPTH = 2000.0 H FREQC = 250.0 HZ 
RNGSTP s 5.04 ANGSTP = 1.0 DEG MRZRNG * 200.0 M YR = 250.0 H 
YERROR * 5.0 H 
EILGENRAYS 
RAY BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETAT(RAD) THE TAR (RAD) TMETAYCRAD) HODULO THIOP 1 


THETAY(RAD) 
l 53.000 251.11 0.167269 52.856 O0.2627462E*03 O.16720491E*03 0.9549710E*02 1.2493187 
2 54.000 245.69 0.165112 53.856 0.20593569E°0S O.169GS32E*03 0.8993371E*02 1.9691119 
3 150.000 246.73 -2£67040 -7. 939 0.4G1°94649E*03 0.1047093E*03 0.3147556E*03 0.5963779 


TOTAL CPU TIME = 0 MIN » 22.47 SEC 


TABLE 3.21 
EIGENRAYS FOUND USING THE 
ODE SOLVER METHOD o 


RAY TRACING USING AKIMA CUBIC SPLINE 8&8 ODE SOLVER 
CASE 2CH4 


NDATA = 4 
ERAYS = T PRTALL © F 
Yo = 100.0 M DEPTH = 2000.0 H FREQC «= 250.0 MZ 
RNGSTP «& 5.04 AHGSTP = 1.0 DEG HRZRHG = 200.0 M YR = 250.0 M 
YERROR = 5.04 
EIGENRAYS 
RAY BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETAT{(RAD) THE TAR CRAD) TMETAYCRAD) MODULO TWOP 


THE TAY(RAD) 
| 53.000 250.85 0.167063 52.050 O.2624GC01K6E*O3 O.16724691E*03 0.9517250E*+02 0.9247239 


2 54.000 245.44 0.164912 Ss oe O.T5904G3SE*O3 O.169G232E*03 =0.8°6S032E 402 1.6557224 
3 150.000 246.56 0.265843 29.971 O.4191563E*03 O.1047093E*03 0.3144470E*03 0.2877732 


TOTAL CPU TIHE = 24 HIN , 12.94 SEC 
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launch angles of 53, 54 and 150 degrees. The tables highlight 
this problem's parameters to be: 

mea SOUFCe depen y, of 100 meters; 

mea, recetver Gepth y of 250 meters; 


- and a horizontal range (HRZRNG) to the receiver of 200 
meters. 


As expected, the ODE solver runs into range limitations 
imposed by CPU time requirements. The 200 meter example 
presented consumed over 24 minutes of CPU time using the ODE 
solver versus the 21 seconds used by the piecewise linear 
method. 

Tables 3.22 and 3.23 (along with their corresponding 
Figures 3.15 and 3.16) show solutions to long-range eigenray 
search problems. These two examples continue to use a source 
Beetney, vet LOO meters and a receiver depth y, of 250 meters. 
Table 3.22 shows a horizontal range (HRZRNG) of 10 kilometers 
watt a depth error y..,. at the receiver of 2 meters. The 
angle step size (angstp) between each ray is 0.1°. Table 3.22 
also shows that three eigenrays were found with launch angles 
of 51.2°, 90.9° and 129.5°. The CPU time for this 10km run is 


shown to be 8.5 minutes. 
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TABLE 3.22 


LONG RANGE EIGENRAYS, 


RAY TRACING USING PIECEWISE LINEAR APPROXIMATIDN : 


CASE 1CH4G 
NDATA = 4G NUMBER OF GRADIENTS = 3 
GC 1) = 0,16000E-01 1/SEC GC 2) = -0.1888°E-01 I1/SEC 
GC 3) = 0.17000E-01 I1/SEC 
ERAYS = T PRTALL = F 
YO = 100.0 ™M DEPTH = 2000.0 M FREOC = 250.0 HZ 
RNGSTP = S.0M ANGSTP = 0.1 DEG HRZRNG 10000.0 M 
YERRDR 3 2.0 M 
ELGENRAYS 
RAY BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETATCRAD) THE TAR (RAD) 
l 51.200 248.91 8.645958 51.067 O.135810G4E+0S 0.8160392E+04 
2 90.900 250.77 6.674996 66.356 0.1048506E*0S 0.1046964E+08 
3 «9.800 249.21 8.730378 $0.370 O0.1371365E+0S5 0.8079626E+04 
TDTAL CPU TIME = 8 MIN + 31.58 SEC 
TABLE 3.23 


LONG RANGE EIGENRAYS, 


RAY TRACING USING PIECEWISE LINEAR APPRDXIMATION 


CASE ICH4 
NDATA = 4G NUMBER DF GRADIENTS = 3 
GC 1) = O.16000E-01 I1/SEC GC 2) = -0.18889E-01 1/SEC 
GC 3) = O.17000E-01 1/SEC 
ERAYS = T PRTALL = F 
YO = 100.0 ™ DEPTH = 2000.0 M FREOC 250.0 HZ 
RNGSTP = $5.0 MM ANGSTPF = 1.0 DEG HRZRNG = 5§0000.0 M 
YERROR = 15.0 M 
E1GENRAYS 
RAY BETAO(DEG) DEPTH(M) TRVLT(SEC) BETAY(DEG) THETATC(RAD) THE TAR (RAD) 
1 40.000 236.69 §2.421823 39.917 0.82546401E*0S O.3365291E+*05 
2 102.000 261.00 34.466009° 100.535 0.5409868E*05 0.8121057E+05 
3 119.000 as0.35 38.520812 60.820 0.608083S5E*0S 0.4579040E+05S 
G4 122.000 257.74 39.734666 122.182 0.6241507E*0S8 0.4439925E+05 
TDTAL CPU TIME = 4 MIN ,» 27.61 SEC 


48 


THE TAY(RAD) 


0.5420648E+04 
0.1842215E*02 
0.5634020E+04 


THE TAYCRAD ) 


0.4869109E*08 
0.2888117E+04 
0.147179SE+08 
0.1801S81E*05 


PIECEWISE LINEAR SOLUTION 


MDDULO TWOP! 


THE TAY (RAD) 
4.$420569 
2.8557634 
4.2855 981 


PIECEWISE LINEAR SOLUTION 


MODULO TWDP1 

THETAY (RAD ) 
2.6894315 
4.13805-¢8 
2.7282732 
1.9203099 
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Table 3.23 shows simulation results using the piecewise 
linear method with a horizontal range (HRZRNG) of 50km. 
Parameters that affect execution time such as angle step sizZe, 
depths, etc., match the simulation presented in Table 3.21 
using the ODE solver. Table 3.21 shows that the ODE solver 
used 24 minutes to complete a 200 meter range eigenvalue 
problem. In comparison, Table 3.23 shows that the piecewise 
linear method completed the 50km run in only 4 minutes. 

In summary, this chapter has shown that 

- the ODE solver performs accurate phase calculations 
with far fewer data pair samples when a_e smooth 
sound-speed profile is encountered; 

- increasing the range step size of the ODE solver to 
five meters reduces the CPU time cost without 


signficantly affecting the phase solution accuracies; 


- the piecewise linear method can run _ long-range 
simulations using relatively little CPU time; and 


- the piecewise linear method requires many sound-speed 


data pair samples for an accurate solution in a medium 
with a smooth, sinusoidal profile. 


Syl 


IV. CONCLUSIONS AND RECOMMENDATIONS 


The transfer function of the inhomogeneous ocean based on 
the WKB approximation requires solving a phase integral. Ray 
acoustics theory can provide solutions to this phase integral. 
The simulations performed show that the travel times 
calculated using the theories of ray acoustics can be used to 
solve the phase integral avoiding direct numerical 
integration. 

Two applications of ray acoustics produced computer 
Simulation codes which 


- are capable of solving for the position, travel time 
and phase of a propagating ray, and 


- have very different advantages and costs. 

The first application was the piecewise linear approximation. 
Sound-speed versus depth data pairs sampled from the ocean 
medium were connected with constant gradient linear segments. 
Well-known, closed form equations form the mathematical model 
for sound propagation. The solutions are low cost (in terms 
of CPU time), but many data samples are required for accurate 
phase solutions for arbitrary sound-speed profiles. 

In contrast, the Akima cubic spline/ODE solver method uses 
the medium samples to form a continuously variable sound-speed 
profile. Accurate phase calculations can be made with a 
minimum of sound-speed versus depth data pairs. The 
disadvantage in using the ODE wee is its exorbitant cost in 
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terms of CPU time. Simulations must be short range problems 
to limit computing costs. 

The ability to search for and identify eigenrays was 
developed for each of the two ray acoustics theory 
applications. Very distinctive characteristics were seen for 
phase calculations in terms of the number of data samples 
needed to assure accurate solutions. The identification of 
eigenrays is an easier task of position or depth computation. 
Here the distinctions tend to disappear with the two methods 
arriving at the same solutions for eigenray launch angles. 
The cost or CPU time required continues to strongly favor 
using the piecewise linear approach. 

These findings indicate that most propagation problems 
will require a piecewise linear approach for computational 
efficiency. Using the ODE solver would severely limit the 
range of computer simulations. 

Careful sampling of the ocean media is required to obtain 
accurate results from the piecewise linear method. Sufficient 
sound-speed versus depth data pairs must be used to accurately 
represent the sound-speed profile. 

This thesis developed the tools to quantify the strengths 
and weaknesses of two phase computation methods in a variety 
of media. With these findings in mind, future work 
recommendations are to 


- incorporate each phase computation technique as a 
module in the larger pulse propagation code; and 


3 


- run pulse propagation simulations to compare the 
received pulse shapes. 


If this comparison shows that the piecewise linear 
approximation produces a relatively undistorted received 
pulse, it is an efficient solution to the phase computation 


problem. 
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